% this is  Si rib waveguide 
% structure is fixed, it is searching modes, next sorting and plot mde
% field distribution

clear all

import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
model.name('AlGaAs_rib_waveguide mode analysis');


SaveToFile='yes';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SaveToFile=='yes'
[MyFile myPath]=uiputfile('*.*','save file','C:\Users\DellBig\Desktop\delete me B');
if MyFile==0;return;end
end

lambda=1.55;

numberOfMode=6;   %  number of meode to look
followPrevious='y';  % use previuos calculations to calculate next

RibWaveGeometry

myChoice=1;    % 1- waveguide width    2 -etching depth  3% core layer

if myChoice==1
ScanningParameter='WireWidth';   
ScanningUnit='[um]';
MyScanningTitle=' waveguide width ';
MyScan=[0.3:0.05:1];

%CalculationWidth=MyScan(length(MyScan))+4.6;
end

           
if myChoice==2
ScanningParameter='EtchingDepth';   
ScanningUnit='[um]';
MyScanningTitle=' EtchingDepth ';
MyScan=0.5:-0.025:0.2;
end

if myChoice==3
ScanningParameter='CoreThickness';   
ScanningUnit='[um]';
MyScanningTitle=' waveguide width ';
MyScan=[0.7 0.65 0.6 0.65 0.5 0.45 0.42:-0.05:0.25];

end





ShowGraph='non';





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.physics.create('emw', 'ElectromagneticWaves', 'geom1');

model.study.create('std1');
model.study('std1').feature.create('mode', 'ModeAnalysis');

model.physics('emw').feature('wee1').set('DisplacementFieldModel', 1, 'RefractiveIndex');
model.physics('emw').feature('wee1').set('n_mat', 1, 'from_mat');

model.study('std1').feature('mode').set('geomselection', 'geom1');
model.study('std1').feature('mode').set('physselection', 'emw');
model.study('std1').feature('mode').set('shift', 'nAprox');  %  aproximate mode refractive index
model.study('std1').feature('mode').set('modeFreq', ['c_const/' num2str(lambda) '[um]']);



model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'mode');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('e1', 'Eigenvalue');
model.sol('sol1').feature('e1').set('shift',  'nAprox');   %  aproximate mode refractive index
model.sol('sol1').feature('e1').set('neigs', numberOfMode);                            %  number of mode to search
model.sol('sol1').feature('e1').set('transform', 'effective_mode_index');
model.sol('sol1').feature('e1').set('control', 'mode');
model.sol('sol1').feature('e1').feature('aDef').set('complexfun', true);
model.sol('sol1').attach('std1');

%model.sol('sol1').runAll;
%disp('calculations OK')
TEFirstPassed='n';TMFirstPassed='n';jf_TM=0;jf_TE=0;


for j1=1:length(MyScan)

model.param.set(ScanningParameter, [num2str(MyScan(j1)) ScanningUnit]);

if myChoice==3  %  make etching depth the same as core
    model.param.set('EtchingDepth', [num2str(MyScan(j1)) '[um]']);
end

if followPrevious=='y'  % make previous solution as initial aprox
if j1>1  %  make previous solution as initial aprox
    if jTM>0
        nAprox=mode(TMmod(1,1)).neff;
    elseif jTE>0
        nAprox=mode(TEmod(1,1)).neff; 
    end;
        model.study('std1').feature('mode').set('shift', num2str(nAprox));  %  aproximate mode refractive index
        model.sol('sol1').feature('e1').set('shift',  num2str(nAprox));
end;
end

        
        

model.sol('sol1').runAll;
MyWidth(j1)=MyScan(j1);
disp(['thickness ' num2str(MyScan(j1))]);

%%% used for integration
MyX=wireWidth/2:SiO2Thickness/20:wireWidth/2+SiO2Thickness;
MyY=zeros(1,length(MyX))+CutLinePos;


% find which solutions are realistic
n_eff_all = mphglobal(model,'emw.neff','Complexout','on');
jMode=0;
for j=1:length(n_eff_all)
    %if real(n_eff_all(j))>nAlGaAs(xClad,lambda)      
        jMode=jMode+1;
        mode(jMode).neff=real(n_eff_all(j));
        mode(jMode).loss=kdb(imag(n_eff_all(j)),lambda)/10;  %loss in dB/mm
       
        mode(jMode).n=n_eff_all(j);
        mode(jMode).solutionN=j;
   % end
end
%disp([' mode total number:  ' num2str(jMode)])

AreaAll=mphint(model,'1','Solnum','1');  % integration over all domains
AreaCore=mphint(model,'1','Solnum','1','selection',model.selection('geom1_sel3').entities(2));
%%%%%% find whether the mode TE or TM
jTM=0;jTE=0;jTEM=0;

for j=1:jMode
    
TEintens=mphint(model,'abs(emw.Ex)^2','Solnum',num2str(mode(j).solutionN))/AreaAll;  % integration over all domains
TMintens=mphint(model,'abs(emw.Ey)^2','Solnum',num2str(mode(j).solutionN))/AreaAll;
mode(j).TEintens=TEintens;mode(j).TMintens=TMintens;
 

if TEintens>3*TMintens 
    TEintensCore=mphint(model,'abs(emw.Ex)^2','Solnum',num2str(mode(j).solutionN),'selection',model.selection('geom1_sel3').entities(2))/AreaCore;
    if 3*TEintensCore>TEintens
        
    jTE=jTE+1;TEmod(jTE,1)=mode(j).solutionN;TEmod(jTE,2)=mode(j).neff;
    end
    
elseif TMintens>3*TEintens
    TMintensCore=mphint(model,'abs(emw.Ey)^2','Solnum',num2str(mode(j).solutionN),'selection',model.selection('geom1_sel3').entities(2))/AreaCore;
       if 4*TMintensCore>TMintens
        
    mode(j).Polar='TM';
    jTM=jTM+1;TMmod(jTM,1)=mode(j).solutionN;TMmod(jTM,2)=mode(j).neff;
       end
else
    mode(j).Polar='TEM';
    jTEM=jTEM+1;TEMmod(jTEM,1)=mode(j).solutionN;TEMmod(jTEM,2)=mode(j).neff;
end

end

% put each mode in order, mode with highest refractive index first
mes='';
if jTE>0 TEmod=sortrows(TEmod,-2); mes=['number of TE modes ' num2str(size(TEmod,1))];end
if jTM>0 TMmod=sortrows(TMmod,-2);mes=[mes ' of TM modes ' num2str(size(TMmod,1))]; end
 if jTEM>0 TEMmod=sortrows(TEMmod,-2);mes=[mes ' of TEM modes ' num2str(size(TEMmod,1))]; end

disp(mes )

if jTE>0
    nTE0(j1)=mode(TEmod(1,1)).neff;
    lossTE0(j1)=mode(TEmod(1,1)).loss;
    
    % find transverse effective refractive index
Ex=mphinterp(model,'Ex','Solnum', num2str(TEmod(1,1)),'coord',[MyX;MyY]);
[kx b2 sigma]=LinearFit(log(abs(Ex)),MyX);
neffxTE0(j1)=abs(kx*lambda/pi/2);
    
else
     nTE0(j1)=NaN;
     lossTE0(j1)=NaN;
     neffxTE0(j1)=NaN;
end

if jTE>1
    nTE1(j1)=mode(TEmod(2,1)).neff;
    lossTE1(j1)=mode(TEmod(2,1)).loss;
else
     nTE1(j1)=NaN;
     lossTE1(j1)=NaN;
end    

if jTM>0
    nTM0(j1)=mode(TMmod(1,1)).neff;
    lossTM0(j1)=mode(TMmod(1,1)).loss;
    
    Ey=mphinterp(model,'Ey','Solnum', num2str(TMmod(1,1)),'coord',[MyX;MyY]);
   [kx b2 sigma]=LinearFit(log(abs(Ey)),MyX);
    neffxTM0(j1)=abs(kx*lambda/pi/2);
    
else
     nTM0(j1)=NaN;
     lossTM0(j1)=NaN;
     neffxTM0(j1)=NaN;
end

if jTM>1
    nTM1(j1)=mode(TMmod(2,1)).neff;
    lossTM1(j1)=mode(TMmod(2,1)).loss;
else
     nTM1(j1)=NaN;
     lossTM1(j1)=NaN;
end  

%%%%%%%%%%%%%%%%%

%%%%%%%%%%%% plot result only zero mode

  if jTE>0
     if TEFirstPassed=='n'
      TEFirstPassed='y';
      dataName='pgTE0';

      model.result.create(dataName, 'PlotGroup2D');
      model.result(dataName).feature.create('surf1', 'Surface');
      model.result(dataName).feature('surf1').set('descr', ['Electric field TE0'  ' mode, x component loss' num2str(mode(TEmod(1,1)).loss) ' dB/mm']);
      model.result(dataName).set('title', [' TE0 mode, abs(Ex) ' ScanningParameter '=' num2str(MyScan(j1)) ' um  n= ' num2str(nTE0(j1)) '  loss: ' num2str(lossTE0(j1)) ' dB/mm']);
      model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ex)');
       model.result(dataName).name(['TE' num2str(j-1)]);
       model.result(dataName).set('solnum', num2str(TEmod(1,1)));
        %model.result('dataName').set('data', 'dset1');
        model.result(dataName).run;
        % model.result(dataName).feature('surf1').set('rangecoloractive', 'on');
        hTE0=figure(1); mphplot(model,dataName,'rangenum',1);
        jf_TE=jf_TE+1;F_TEO(jf_TE)=getframe(hTE0);

     else
              
      dataName='pgTE0';
      model.result(dataName).set('solnum', num2str(TEmod(1,1)));
      model.result(dataName).set('title', [' TE0 mode, abs(Ex) ' ScanningParameter '=' num2str(MyScan(j1)) ' um  n= ' num2str(nTE0(j1)) '  loss: ' num2str(lossTE0(j1)) ' dB/mm']);
      model.result(dataName).run;

        hTE0=figure(1); mphplot(model,dataName,'rangenum',1);
        jf_TE=jf_TE+1;F_TEO(jf_TE)=getframe(hTE0);
      
     end 
end

if jTM>0
     if TMFirstPassed=='n'
      TMFirstPassed='y';
      dataName='pgTM0';
model.result.create(dataName, 'PlotGroup2D');
model.result(dataName).feature.create('surf1', 'Surface');
model.result(dataName).feature('surf1').set('descr', ['Electric field TE' num2str(j-1) ' mode, x component']);
model.result(dataName).set('title', [' TM0 mode, abs(Ey) ' ScanningParameter '=' num2str(MyScan(j1)) ' um n= ' num2str(nTM0(j1)) '  loss: ' num2str(lossTM0(j1)) ' dB/mm']);

model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ey)');
model.result(dataName).name(['TM' num2str(j-1)]);
model.result(dataName).set('solnum', num2str(TMmod(1,1)));
%model.result(dataName).set('data', 'dset1');
model.result(dataName).run;
%model.result(dataName).feature('surf1').set('rangecoloractive', 'on');
        hTM0=figure(2); mphplot(model,dataName,'rangenum',1);
        jf_TM=jf_TM+1;F_TMO(jf_TM)=getframe(hTM0);
     else
         dataName='pgTM0';
         model.result(dataName).set('title', [' TM0 mode, abs(Ey) ' ScanningParameter '=' num2str(MyScan(j1)) ' um n= ' num2str(nTM0(j1)) '  loss: ' num2str(lossTM0(j1)) ' dB/mm']);
         model.result(dataName).set('solnum', num2str(TMmod(1,1)));
         model.result(dataName).run;
         hTM0=figure(2); mphplot(model,dataName,'rangenum',1);
        jf_TM=jf_TM+1;F_TMO(jf_TM)=getframe(hTM0);
         
       end 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (jTM>0)&&(jTE>0)
    PhaseShift0(j1)=(mode(TMmod(1,1)).neff-mode(TEmod(1,1)).neff)*360/lambda;  % deg/um
else
    PhaseShift0(j1)=NaN;
end

if (jTM>1)&&(jTE>1)
    PhaseShift1(j1)=(mode(TMmod(2,1)).neff-mode(TEmod(2,1)).neff)*360/lambda;  % deg/um
else
    PhaseShift1(j1)=0;
end

 





end


if SaveToFile=='yes'

   
   
 fid = fopen([myPath MyFile '.txt' ],'w');
 
    fprintf(fid,[ScanningParameter ', TE0_n, TE0_loss, TE0_nx, TM0_n, TM0_loss, TM0_nx, PhaseShift0, TE1_n, TE1_loss, TM1_n, TM1_loss, PhaseShift1\n']);
    %fprintf(fid,['um no  dB/mm no  dB/mm deg/um no  dB/mm no  dB/mm deg/um  \n']);
    
    for j2=1:length(MyScan)  
        MyPrint=[num2str(MyScan(j2)) ];
        tmp=[nTE0(j2) lossTE0(j2) neffxTE0(j2) nTM0(j2) lossTM0(j2) neffxTM0(j2) PhaseShift0(j2)];
        tmp=[tmp nTE1(j2) lossTE1(j2) nTM1(j2) lossTM1(j2) PhaseShift1(j2)];
        for j=1:length(tmp)
            MyPrint=[MyPrint ','];
            if ~isnan(tmp(j))
                MyPrint=[MyPrint num2str(tmp(j))];
            end
        end

       fprintf(fid,[MyPrint '\n']);
    end
    
 fclose(fid); 
 
 if TEFirstPassed=='y'
     F1=F_TEO;
     save([myPath MyFile 'TEO'], 'F1');  %save images frame
 end
 if TMFirstPassed=='y'
     clear F1; F1=F_TMO;
     save([myPath MyFile 'TMO'], 'F1');  %save images frame
 end


 
end

h4=figure(4);movie(h4,F_TMO,100,1);

disp('All OK')